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SUMMARY 


Approximate nonlinear inviscid theoretical techniques for predicting 
aerodynamic characteristics and surface pressures for relatively slender 
vehicles at supersonic and moderate hypersonic speeds were developed. 
Emphasis was placed on approaches that would be responsive to conceptual 
configuration design level of effort. Second order small disturbance and 
full potential theory was utilized to meet this objective. 

Numerical codes were developed for relatively general three 
dimensional geometries to evaluate the capability of the approximate 
equations of motion considered. Results from the computations indicate 
good agreement with experimental results for a variety of wing, body, and 
wing-body shapes. 
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1. INTRODUCTION 

An examination of the literature for supersonic/hypersonic aircraft 
provides an indication of the flexibility and generality required for a 
prediction technique. Typical configuration development variables include 
wing section, incidence, height, dihedral, planform, effectiveness of 
longitudinal control surfaces for trim, effectiveness of empennage for 
directional stability, and propulsion system-airframe interactions. 

State-of-the-art response to these prediction requirements is provided 
by hypersonic impact methods as well as linearized analysis and design 
algorithms. These approaches can treat complex geometries with minimum 
response time and cost, with efficient predicted data coverage in terms 
of Mach number, angle of attack, trim deflection, yaw angle, etc. Shortcomings 
are present, however, in both the impact and linearized methods. For the 
former, interference between surface elements is totally ignored in 
implementations such as classicial Newtonian, tangent wedge, and cone theories. 
Cross- flow interactions and stagnation point singularities are also implicitly 
disregarded. In the latter, shocks, vorticity, and entropy wakes and layers 
are excluded. Furthermore, superposition of elementary solutions such as 
those for thickness and angle of attack freely used in linear models are, 
strictly speaking, invalid at hypersonic speeds. 

A need exists for new aerodynamic prediction techniques to optimize 
vehicles designed to travel at supersonic/hypersonic speeds. One requirement 
of a new aerodynamic prediction technique is that it be more accurate than 
simple noninterfering panel methods. Another specification is that it be more 
computationally efficient than currently available explicit finite -difference 
methods so that it can be incoporated into a practical design procedure. The 
new approach should include enough of the physics of the flow to allow realistic 
optimization and should permit consideration of appropriate interactions 
between components of promising arrangements, since this has been found to be 
the key to increasing aerodynamic efficiency using linear methodology. 

Nonlinear potential theoretical formulations hold the promise of meeting this 
objective and providing economic design codes which are responsive to 
conceptual vehicle definition efforts. 
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2. LIST OF SYMBOLS 


AR 

b 

c 

c 

CD 

Cl 

Cm 

C P 

c z 

d 

D 

F 

l 

L 

M 

P 

q 

s 

u,v,w 

x,y,z 

a 

6 

e 

<P 

[ 3 
x 

A 


aspect ratio, b^/S 

wing span 

chord 

mean aerodynamic chord 
drag coefficient, D/q^S 
lift coefficient, L/q^S 
pitching moment coefficient, M/q^Sc 
pressure coefficient, P-Poo 
normal force coefficient, Fy/q^S 
body maximum diameter 


drag 

force 

body length 
lift 

Mach number or moment 
static pressure 
dynamic pressure 
reference area 

axial, vertical, and lateral perturbation velocity 

axial, vertical, lateral body axis cartesian coordinate 

angle of attack 

sideslip angle 

see figures 5.2 and 5.3 

polar angle 

velocity potential 

denotes jump 


taper ratio, 



sweep 
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SUBSCRIPTS 

ac 

LE 

r 

t 

x,y,z 

a 


aerodynamic center 
leading edge 
root 

tip or maximum thickness 

component in x,y,z direction respectively 

angle of attack 

frees tream 
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3. METHODOLOGY 


Emphasis is placed on approximate theoretical approaches which are capable 
of treating relatively general three dimensional problems but still sufficiently 
simple to be responsive to vehicle conceptual design efforts . The basic intent 
of the methodology is to produce improvement in lift-drag ratio of supersonic/ 
hypersonic cruise vehicles. As a result of the strong impact that favorable 
interference has had on supersonic design and the use of such concepts in recent 
advanced hypersonic aircraft studies, candidate analysis should be general 
enough to systematically treat such problems. Finally, interest in high 
aerodynamic efficiency emphasizes relatively slender configurations at modest 
attitude; that is, moderate values of the hypersonic similarity parameter. 

Prior theoretical effort ^ ^ has advanced the supersonic/hypersonic 
aerodynamic prediction state-of-the-art at the conceptual design level. 
Numerical second-order potential small disturbance analysis was developed as a 
first step up from the widely used linear theory. Such a formulation 
incorporates nonlinear behavior by iteration of the Prandtl-Glauert 
approximation. This approach is known to extend the prediction success for 
airfoil and cone surface pressure to substantially higher values of the 
hypersonic similarity parameter than the first-order theory. The next level of 
theoretical richness vis-a-vis a full -potential equation of motion formulation 
was explored in parallel. This analysis eliminates edge singularities and 
improves treatment of characteristic surfaces but still retains isentropic 
assumption. 

Hypersonic small disturbance theory was considered in an earlier study! 
in recognition of the progressive non- isentropic behavior of the flow as the 
value of the hypersonic similarity parameter increases. Finite difference 
analysis of this approximation 4 indicated that the solution was essentially 
as complex as that for the Euler equation and thus would not be particularly 
responsive to conceptual design level of effort. This approach was not 
pursued in the present study on the basis of this finding and the previously 
cited success of potential analysis at moderate hypersonic conditions. 
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4. SECOND ORDER DESIGN/ANALYSIS 
1 2 3 

The second-order analysis ’ and design codes developed under this contract 
are based on correcting the three-dimensional first-order potential solution, <J>, 
for the dominant contributions of the nonhomogeneous Prandtl-Glauert equation 




[( i - mL . ) * . (« I *..) 2 


where 
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32 


(1-MJ + + \ * 
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and associated surface boundary conditions. The local surface pressure 
coefficient has the fonn 


c - c ^ 

p p p 


where the superscripts refer to the first- and second-order contributions, 
respectively. 

4.1 GEOMETRY 

The Aerodynamic Preliminary Analysis System (APAS)^*^ was interfaced with the 
second order analysis code to achieve general configuration modeling since it 
provides a well automated interactive capability. Modifications of variables 
such as surface planform, component location, incidence, dihedral, etc are easily 
accomplished through use of system interactive commands. A typical APAS 
aircraft model is presented on figure 4.1 to illustrate the geometric generality 
available. 


4.2 INPUT/OUTPUT GRAPHICS 


The responsiveness of aerodynamic analysis at the conceptual design level 
is dependent on the ability to verify modeling and review/ col late solution 
results for a large number of parametric configuration variables jm a timely 
manner. The use of graphics to display and examine geometry paneling 
minimizes modeling errors and associated erroneous analysis. Typical graphical 
finite element configuration definition is presented on figure 4.2 for a wing- 
body- canard arrangement. 

Review and comparison of t^ie local characteristics is used to assess the 
impact of parametric geometry modification. Chordwise distribution of surface 
pressure and the spanwise variation of typical sectional characteristics such 
as lift, drag, leading edge thrust, center of pressure, etc may be displayed 
for this purpose. A typical set of output is presented on figure 4.3. 
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A PAS 


^rcraft Model 
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a) Plan View 



b) Orthographic View 


Figure 4.2 Typical Finite Element Model Graphics 
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a) Chordwise Pressure Distribution 
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b) Sectional Lift 


[ypical Aerodynamic Analysis Graphics 
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d) Sectional Center of Pressure 
F igure 4.3. Cone luded 
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5. FULL POTENTIAL ANALYSIS 


The full potential equation in conservative or nonconservative form is 
frequently used for transonic flow analysis, where the local Mach number does 
not exceed approximately 1.4. If the assumptions of irrotationality and 
isentropicity are reasonably valid, then the full potential equation is 
expected to yield results comparable to Euler equations, even for supersonic/ 
hypersonic flow fields. For conceptual design studies, where short response 
time is desired, the full potential methods can be an attractive substitute for 
expensive Euler methods and less accurate linear theory methods. 

A nonlinear aerodynamic prediction technique based on the full potential 
equation in conservation form has been developed, during this contractual 
effort, for the treatment of supersonic flows. A detailed description of the 
method has been presented in several published papers. The most recent are 
enclosed in Appendix A for convenience. The first paper describes the method 
for the treatment of predominantly supersonic flows with regions of subsonic 
flow that usually occur at low supersonic Mach numbers The second paper 
describes the grid generation method that is used in the full potential code 
to treat general aircraft configurations. Since the Appendix describes the 
full potential method, only the material not included in the published articles 
will subsequently be presented. 

5.1 GEOMETRY 


Configuration geometry input to the full potential code is described in 
reference 8 and summarized briefly here. The geometry is prescribed as a set 
of discrete points at various axial cross-sectional locations. These input points 
are usually obtained from a geometry package such as GEMPAK^ or CDS^O . The input 
points are then divided into several groups or patches, and for each patch a 
key-point system is established. The geometry at any marching plane is then 
obtained by joining the appropriate key-points for each patch as shown in figure 
5.1. A cubic spline fit of the key-points is used to obtain the grid point 
distribution (clustering) on the body surface. The physical grid between the 
body surface and the outer free stream boundary is then generated using an 
elliptic grid solver. 

5.2 WAKE MODEL 

The first paper of the Appendix reported the development of a wake treatment 
for the full potential method. This treatment has now been extended for the case 
of swept trailing edges. Figure 5-2 shows the schematic of a general swept 
trailing edge wake system. As stated in the paper, in order to treat the region 
behind the trailing edge, an artificial cut is created and the pressure jump, [P] , 
across such a cut is imposed to be zero as a boundary condition. (Note the 
nomenclature is that used in the Appendix) . This is achieved by 
maintaining the jump in the velocity potential, <p, along a K = constant line (see 
figure 5.2) for J => 2 to be the same as the jump in <j>, [<f>] , at the trailing edge. 
The full potential equation is not solved at grid points on the wake. cut. Instead 
the second derivative of <p with respect to the radial direction n, 4>nn> equals 
zero is solved to provide [<J>n] = 0 across the wake cut. Maintaining [<j>] constant 
along a K line provides [<jvj = 0. The combination of = 0 and [cj>ri] = 0 
across the cut approximately satisfies [P] =0. 
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5.3 SIDESLIP 
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A complete analysis of an aircraft configuration must address the vehicle's 
performance under angle of attack (a) and sideslip (6) flight conditions. For 
asymmetric configurations a similar pattern exists in that the flow field is not 
symmetric and therefore required a solution of the entire crossflow plane. 
Reference 11 describes the addition of sideslip to the full potential code and 
the highlights are given below. 

For pitch problems, only the half plane needs to be solved with the plane 
of symmetry boundary conditions imposed along K = 2 and (KMAX-1) , as shown 
in Figure 5.3a. Imposing the flow conditions along K = 1 are to be the same 
as the ones along K = 3, the Lr operator results in a tridiagonal system that 
can be easily solved. When sideslip is present the entire cross-flow plane 
needs to be solved as shown in Figure 5.3b. In this case, the flow conditions 
along K = 1 are set to be the same as the ones along K = KMAX-2. This destroys 
the tridiagonal nature of the Lg operator. A special routine has been 
developed to invert a matrix of the following type. 


— X 0 

X X 'x — o o 

' \ S 

* X X 

Vx ' 

' ' V 

\ v N 

Co X X X 

OX 


CD 


In the current formulation, positive angle of attack a represents a 
positive Cartesian velocity v in the free stream and similarly positive 
sideslip 6 produces a positive w in the free stream. When both angle of 
attack and yaw are present , first the free stream is turned by an angle 
6 and then by a. Let U,y,z) be the inertial Cartesian system. After 
an initial sideslip turn 8 let the wind axis^system be (x’,y',z'), and 
after a subsequent a turn let it become (x,y,z). Then, reterring to 
Figure 5.4 
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li 



The free stream is now along x. The normalized free-stream velocity potential, 
O is given by 


= 4>oo = x cos a cos 0 + y sin a + z cos a sin 0 (3) 

Uco 

Using equation (2), the drag, lift and side force are 

Drag = F x cos a cos 0 + F y sin a + F z cos a sin 0 

( 4 ) 

Lift = F x sin a cos 0 + Fy cos a 
Side Force = -F x sin 0 + F z cos 0 
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Figure 5.1. Full Potential Geometry Input 
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Figure 5.2 Swept Trailing Edge Wake Analysis 
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K = 1 2 3 



a) PLANE OF SYMMETRY 

(8 = 0 ) 

<t>l 1 = 3 

^j.KMAX * <t>\, KMAX-2 


KMAX-1 



b) PERIODIC CONDITION FOR 

(8 * 0 ) 

= ^j, KMAX-2 
<*>j, KMAX = ^j,3 


Figure 5.3. Boundary Condition Treatment 
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Figure 5.4. Notation for Sideslip and Angle of Attack 
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6.0 RESULTS 


6.1 SECOND ORDER 


A slender delta wing-body test case? was initially considered for code 
validation purposes. The general arrangement is presented on figure 6.1.1. 

The wing had a leading edge sweep of seventy degrees and a maximum thickness of 
four percent based on sectional chord length. The body fineness ratio, £/d, was 
twelve with a body width to wing span ratio, w/b, of twenty five percent. 
Experimental data and first and second order predictions covering a Mach 
number range of 2 to 7.4 are presented in figure 6.1.2 for lift curve slope, 
longitudinal stability parameter, dCm/dCL, and drag due to lift factor*. The 
experimental test results were generally linear up through maximum lift-drag 
ratio and consequently the previously cited aerodynamic parameters are constant 
in agreement with second order theory. Comparison of the estimations with 
measurements indicate a systematically improved correlation for second order 
analysis. The general prediction success is considered to be satisfactory 
for conceptual /preliminary design studies. Further, the analysis is responsive 
to such activities since a total CPU time of two minutes on an IBM 3033 was 
required to generate the results for the seven Mach numbers considered. 


*Minimum drag was not correlated as a result of uncertainties associated with 
prediction of untripped experimental friction drag levels. 
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Figure 6.1.1. Delta Wing-Body Geometry 
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Figure 6,1.2. Delta Wing-Body Test Case 
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6.2 FULL POTENTIAL 


Figure 6.2.1 presents the pressure distribution on the surface of an arrow- 
wing configuration 1 ' 2 at axial location x/SL = 0.8 for = 2.98 and a = 10.01°. 
The improvement in the prediction using the wake model is illustrated. The 
dashed line represents the result for no wake treatment (assumes a flat plate 
behind the trailing edge) and the solid line shows the improvement in the pressure 
distribution for a zero jump in pressure across the wake cut. The solid line 
pressures on the body agree very well with experiments and are expected to 
improve force and moment calculations. 

Figure 6.2.2 presents the pressure distribution on a forebody 1 ^ at M m = 2.5 
and a sideslip 8 = 5.02°. The axial pressure distributions compare well with 
experimental data. Figure 6.2.3 shows the circumferential pressure distribution 
for the same forebody at x/l = 0.68 for = 1.7, a = 10° and 6 = 5.02°. The 
comparison with experimental data is again very good. 
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Figure 6.2.1. Pressure Distribution at x/i = 0.8 for an Arrow 
Wing; M ot = 2.96, a = 10.01° 
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X/l 


Figure 6.2.2. Pressure Distribution on a Forebody in Sideslip 
M = 2.5, 6 = 5.02°, a = 0° 
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Full Potential 



Figure 6. 2.3. Circumferential Pressure Distribution at x/l = 0,68 
for a Forebody in Pitch and Sideslip, Mo, = 1.7, 
a = 10.0°, 6 = 5.02° 
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7. CONCLUSIONS 


Based on the theoretical development and comparison with experimental 
measurements described in this document, the following conclusions are 
made. 

1. Improved prediction of supersonic/hypersonic aerodynamic 
characteristics and surface pressures for general wing-body 
shapes has been demonst rated using nonlinear potential 
analysis for values of the hypersonic similarity parameter, 

M6, less than one. 

2. Full potential analysis successfully eliminates subsonic edge 
singularities and linear characteristics approximations of 
second order theory. The formulation is an order of magnitude 
faster than Euler solvers while maintaining comparable 
prediction accuracy for flows of moderate shock strength. 

3 

Nonlinear potential theory provides advanced aerodynamic 
prediction techniques that are responsive to the 
conceptual design problem at supersonic and moderate 
hypersonic conditions. 
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Abstract 

A nonlinear method based on the full potential 
equation In conservation form, cast in an arbitrary 
coordinate system, has been developed to treat pre- 
dominantly supersonic flows with embedded subsonic 
regions. This type of flow field occurs frequently 
near the fuselage-canopy junction area and wing 
leading edge regions for a moderately swept fighter 
configuration. The method uses the theory of char- 
acteristics to accurately monitor the type- 
dependent flow field. A conservative switching 
scheme is developed to transition from the super- 
sonic marching algorithm to a subsonic relaxation 
procedure, and vice versa. An implicit approximate 
factorization scheme is employed to solve the 
finite-differenced equation. Results are shown for 
a few configurations, including a wing-body-wake 
realistic fighter model having embedded subsonic 
regions. 


I. Introduction 

Nonlinear aerodynamic prediction methods based 
on the full potential equation are used regularly 
for treating transonicl»2 and supersonic 3-3 flows 
over realistic wing-body configurations. The 
transonic algorithms ^ »2 are designed to treat pre- 
dominantly subsonic flows with pockets of super- 
sonic regions bounded by sonic lines and shocks. 

The supersonic methods 3-3 are based on a marching 
concept, and require the flow to remain supersonic 
in a given marching direction. Once the marching 
direction velocity becomes subsonic, the domain of 
dependence changes, and a pure marching scheme 3-3 
will violate the rules of characteristic signal 
propagation. The possibility of a marching veloc- 
ity becoming subsonic in a supersonic flow is 
great, especially for low supersonic freestream 
Mach number flows (M^ = 1.3^ 1.7) over moderately 
swept fighter-like configurations (sweep angle 
A = 45n,50°), and over forebody shapes having a 
sizeable fuselage-canopy junction region. There is 
a strong need to construct a supersonic marching 
computer program that has built-in logics to detect 
and treat the embedded subsonic regions. 


‘Manager, Computational Fluid Dynamics Group, 
Associate Fellow AIAA 

“Member Technical Staff, Member AIAA 

+ Professor, Department of Mathematics, Member AIAA 


The method of Ref. 5 is based on the character- 
istic theory of signal propagation and uses a gen- 
eralized, nonorthogonal , curvilinear coordinate 
system. Compared to other nonlinear supersonic 
methods 3 , the method of Ref. 5 has no restrictions 
(limitations of the full potential theory hold) on 
its applicability to complex geometries and intri- 
cate shocked flow fields. It is a conservative 
formulation and uses numerical mapping techniques 
to generate the body-fitted system. The purpose 
of this paper is to describe an extension to the 
methodology of Ref. 5 to include the treatment of 
embedded subsonic regions in a supersonic flow. 

The paper describes the characteristic theory 
involved in determining the condition for a march- 
ing direction to exist. Once that condition is 
violated, the marching scheme is transitioned to a 
relaxation scheme through a conservative switching 
operator. For marching condition violation, the 
total velocity q does not have to be subsonic. 

Even for a supersonic total velocity q, if the 
component in the marching direction is subsonic, 
a relaxation scheme is required. In order to 
properly produce the necessary artificial viscosity 
through 1 density biasing, the paper defines two 
situations: one, the total velocity q is super- 

sonic, but the marching direction component is sub- 
sonic (defined as Marching Subsonic Region (MSR) in 
the paper), and two, the total velocity q is sub- 
sonic (termed as Total Subsonic Region (TSR) in 
the paper) . 

Results are presented for a few configurations 
that exhibit either the MSR or both the MSR and 
TSR flow field. The paper also presents results 
from a wake model applied to a realistic wing-body 
fighter configuration. 

The Appendix describes a flux biasing concept 
which will supersede the density biasing procedures 
currently in use. 

The methodology of this paper is not restricted 
to the full potential equation alone. Currently, 
similar marching/relaxation methods are under 
development at Rockwell for application in parabol- 
ized Navier-Stokes (PNS) codes to treat the 
embedded subsonic regions or streamwise separated 
flows without having to use a time-dependent 
Navier-Stokes program. 


II. Equation and Characteristic Theory 

The conservative full potential equation cast 
in an arbitrary coordinate system defined bv 
C = c(x,y,z), n = n(x,y,z) , and £ = E(x,y ,z) , 
takes the form 
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where U, V, and W are the 
components. Introducing 
for convenience 

contravariant 
the following 

velocity 

notation 

u = u 1 

, v = u 2 , 

w = u 3 


X * X.| 

. y = * 2 , 

z = x 3 


X 

II 

AJ* 

» H 8 X2 * 

E = X 3 



the contravariant velocities and density are 
given by 



i = 1.2,3 

i = 1,2,3 (transforma- 
j = 1,2,3 tion metrics) (2) 


\/c ( ''" i) /m£ . 


The Jacobian of the transformation J is represented 
by 


, _ atc.n.O - 

J '3&- 




«U 


(3) 


r-j (pw V 


j <*), 


1 


0 

-1 



The subscripts in Eq. (4) denote differentiation 
with respect to that variable. 

The matrices A, B, and C appearing in Eq. (4) 
can now be analyzed to determine the character of 
that equation. In general, the following is true: 

1) Equation (4) is elliptic in the £ direction 

if the matrix A'MaB+BC) has complex 
eigenvalues for all combinations of a and B 
such that a 2 + B 2 = 1 . 

2) Equation (4) is hyperbolic in the £ direc- 
tion if A _1 (aB+BC) has real eigenvalues 
for all a and 8 satisfying a 2 + B 2 = 1. 

The eigenvalue structure of A'MclB + BC) can be 
Obtained by setting the determinant 

| ctB + BC - XA | = 0 (assuming A' 1 exists) . (5) 

Substituting for A, B, and C from Eq. (4), the 
eigenvalues of Eq. (5) are given by solving the 
quadratic 


Equation (1) is in terms of a general coordinate 
system (t.n.E) and can acconmodate any kind of 
mapping procedure, either analytical (conformal 
mapping) or numerical type. The nature of Eq. (1) 
can be analyzed by studying the eigenvalue system 
of Eq. (1). Combining the irrotationality condi- 
tion in the (£,n) and (£,5) plane and Eq. (1), one 
can write the following matrix equation 


\ a( 


-X 2 (pU) A + X|a(pV) + B ( pW ) - + a(pU) + a(pU) 
£ 


- |a 2 (pV), + B 2 (pW). + a8 (pWL + (pV) 

4>_ M>r <P_ H 


( 6 ) 
= 0 . 


Representing Eq. (6) in the form 


Af + Bf + Cf r = 0 
C n £ 


(4) 


AX 2 + BX + C = 0 , 


(7) 


where 

1 0 

0 1 

J <»*>♦ j { °x 1 

0 0 

0 0 - 


the discriminant (B 2 - 4AC) determines the character 
Of Eq. (4): 

1) If (B 2 - 4AC) remains positive for all o and 
B satisfying a 2 + B 2 = 1, then the eigenval- 
ues of Eq. (4) are real and direction £ is 
hyperbolic (marching scheme is valid). 

2) If (B 2 -4AC) is negative , then the eigen- 
values of Eq. (4) are complex and direction 
£ is elliptic (requires a relaxation 
method) . 

To analyze when the eigenvalue solutions of Eq. (6) 
are real and when complex, the discriminant 
(S 2 - 4AC) is rewritten in the following form 
(using Eq. (2)) 




B = 


-1 
0 
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Using the properties of a positive definite quadra- 
tic form and the Schwarz inequality t a -j -j a j j > a ij)« 

Eq. (8) can be shown to have the following results: 

1) (B 2 -4AC) is positive if is less 

than zero. Then the c direction is hyper- 
bolic (the marching algorithm of Ref. 5 
is val id) . 

2) ( B 2 - 4AC) is negative if is 

greater than zero. Then the c direction 
is elliptic (requires a relaxation scheme). 


A. Physical Interpretation 


The physical interpretation of these results 
from the characteristic theory is illustrated in 
Fig. 1. Let q be the total velocity. The projec- 
tion of q in the direction normal to the c=constant 
surface is given by 


q*n = (ui + vj + w(c) 


c 2 + c 2 + C 2 

x y z 


(9) 


v a ll 

where u, v, and w are the Cartesian velocities and 
n is the normal to the c=constant plane. Figure la 
shows the case when — ~ ■ is greater than the speed 


of sound 



For this case, the char- 


acteristic cone of influence is behind the 
C=constant plane and marching along c is valid. 
Figure lb illustrates the case for the q > a, but 




< a 


situation For 


case, a part of the characteristic cone of influ- 
ence lies forward of the c=constant plane and 
marching along £ is not possible. This case 
(Fig. lb) is termed Marching Subsonic Region (MSR) 
in this paper. Figure 1c shows the case when q < a 



This represents a 


pure subsonic flow and marching along c is not pos- 
sible. This case is termed Total Subsonic Region 
(TSR). For cases represented in Figs, lb and lc, 
a relaxation algorithm is required. 
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q ■= TOTAL VELOCITY 
n » NORMAL TO 

t * CONSTANT PLANE '-CHARACTERISTIC 

8 - SPEED OF SOUND ZONE 0F INFLUENCE 

U 2 

b) MARCHING SUBSONIC REGION (MSR). <a 11 - — ) >0.q > a. 
MARCHING ALONG t IS NOT VALID. ® 2 



Fig. 1. Role of characteristics in defining 
supersonic region, Marching Subsonic 
Region (MSR), and Total Subsonic 
Region (TSR). 


III. Numerical Method 


Figure 2 shows the schematic of a fuselage- 
canopy forebody geometry with an embedded MSR and 
TSR present in a supersonic flow. To solve this 
problem, the marching scheme of Ref. 5 will be used 


when 




UPSTREAM 
COMPUTATIONAL 
BOUNDARY FOR 
THE RELAXATION 
ZONE 


s negative, and a relaxation scheme 

SC83 -22339 

u 2 

q> a; (a^ - — ) >0 
c (MARCHING SUBSONIC REGION) 

DOWNSTREAM 
BOUNDARY FOR 
THE RELAXATION 
ZONE 



(TqiAL_SUBSONIC_R EG ION )_ 
MARCHING REGION- 


Fig. 2. Embedded subsonic bubble in a super- 
sonic flow. 
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when 


(■11 -S) 


is positive. First, march from the 


nose up to the plane denoted by (A-B) in Fig. 2, 
ng the method of Ref. 5. Then, between (A-E 
and (C-D), which embed the subsonic bubble (MSR 


using the method of Ref. 5. Then, between (A-B) 


and TSR), use a relaxation scheme and iterate until 
the subsonic bubble is fully captured. Then, 
resume the marching scheme from the plane (C-D), 
downstream of the body. 

The purpose of this paper is to present a con- 
servative algorithm that will automatically switch 
from a pure marching scheme of Ref. 5 to a relaxa- 
tion method at the onset of an MSR formation and 
revert to the marching procedure when the flow 
becomes fully supersonic again. The entire flow 
field can be classified into three types with 
respect to the marching direction C: 

1) At a grid point, the marching direction is 
hyperbolic and the total velocity q is 

supersonic, < 0, q > a. This 

point will use the algorithm of Ref. 5. 

2) At a grid point, the marching direction ; 
is elliptic , > 0* the total 

velocity q is supersonic, q > a (MSR). 

This point will be treated by a transonic 
operator with a built-in density biasing 

based on the magnitude of 

3) At a grid point, the direction c is 
elliptic and the total velocity q is 
subsonic, q < a (TSR). This point will 

be treated by a subsonic central differenced 
operator. 

A. Treatment of |p in Eq. (1) 

Refer to the computational molecule in Fig. 3. 


a 

3c 


B) 



*i 3c 


supersonic 

+ O- 6 ^) 



(10) 


marching 

subsonic 

where 

$ refers to backward differencing 
$ refers to forward differencing 
1 if la,, -^-| < 0 


6. 

1 


= 0 i 


(*.i-5) 
,f ( a " ' 7) 


> 0 . 


In Eq. (10), the first term corresponds to the 
supersonic marching operator of Ref. 5, and the 


1<0 


i ♦ 1 
j 


i-2 i - 1 i ■ * 1 i ♦ 2 i ♦ 3 

fi-i 


6 : * 1; 5; 4. « ■ 1 


6T'f>-T; ( ‘jVi IFROMEQ.10I 

9) SUPERSONIC UPWIND 
OPERATOR 


MARCHING 
SONIC 
LINE — v 

U 2 V U 2 
<• 11 — 2 > <0 /< B 11 ~J 2 » >0 

8 / t i + 1 

— u i T 
— . — +J— 
i-2 i - 1 1' i+ 1 i + 2 


fi-i 
♦ 1*0 

i. i«yv. i ,pjj i 

4i r J 1 n 1 j >* * i * ST 1 r'i + i 


*i-i- 

A 1 _ A AM i , 6 .pV 


bJ BEGINING OF SUBSONIC 
BUBBLE. 


MARCHING 

SONIC 

U 2 y-UNE 

« J \ ‘*11 — 2* < 0 
U t \ _ J2-. 

i-1 i i ♦ 1 i ♦ 2 i + 3 
\ 

I 


— t°— ) » 0 
5[ 1 J ' 

c) END OF SUBSONIC BUBBLE. 


U 2 

(.1,— ,l>0 


i-1 i I + 1 i + 2 i + 3 
*i*i; »i + 1 ”0 

J ±,e. U, 
n'j 1 n'j >i*i 

dl PURE SUBSONIC REGION 

CENTRAL DIFFERENCE OPERATOR 
(c ( IS BACKWARD DIFFERENCEDI 


Fig. 3. Conservative type-dependent switching 
scheme for the treatment of subsonic 
bubble in a supersonic flow. 

second term is the subsonic operator. 

The backward difference operator in Eq. (10) is 
represented by 

w ( p $),„ 

+ /a . yv\ i_ 

+ / a _W\J_ 

\ 13 >2 /« 

A<t> = U i+ i • 


A<(> 


Aif> + U. 


(ID 


The term •— (Atf>) is backward differenced. 

Reference 5 gives more details on this supersonic 
marching operator. 

The forward difference operator in Eq. (10) is 
represented by 


f I 

g 


c +a l2*n 


+ a 13^ 


where 


<rv.-^? + i -^* forU>0 


v = max 




( 12 ) 


(13) 
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The superscript n+1 denotes the current relaxation 
cycle for a subsonic bubble calculation. 


Note that in Eq. (12) the term <4 is backward 

’ 

differenced such that ^ will provide the 


central differencing needed for an elliptic 
(subsonic) point. The density biasing, Eq. (13), 
is activated only when the total velocity q is 
greater than the speed of sound a. This will take 
place when a grid point is in the region denoted by 
MSR in Fig. 2. When q < a, region TSR in Fig. 2, 
the density is not biased and the generation of 
artificial viscosity is turned off. The <f> deriva- 
tives in Eq. (l3) can be rewritten in terms of 
just like in Eq. (11). 
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e 


i+l 


= 0 



(MSR) . 


When 6 i+1 = 1, that is, the point is supersonic 

with respect to c, only the first term in Eq. (15) 
is used and the biased density e is defined by 
(for V > 0) 

^j+H = {1 " ~j +! s >p j + H + 1 W p j + p j-l } (16) 
where 

v = max^O, l-a 22 . 


Equation (10) can also be interpreted as 



elliptic 

operator 


- Vi | (° • 04) 

flux biasing to produce 
the artificial viscosity 


In Eq. (16), the evaluation of p* depends on 
whether the flow is conical or nonconical. For 
conical flows, all p* quantities are evaluated at 

the i th plane. For nonconical flows, at each non- 
conical marching plane, initially p* is set to be 

the value at the i th plane and then subsequently 
iterated to convergence by setting p* to the previ 
ous iterated value of p at the current i+l plane. 
Reference 5 provides more details on the density 
biasing procedure and the implicit treatment of 

l (' J) M E <- " 5) - 


Figure 3 illustrates various possibilities that 
can be handled by Eq. (10). It has both the shock 
point operator and the sonic operator required to 
treat the type-dependent flow. The only issue that 
philosophically affects the concept of a conserva- 
tive scheme is that the definition of pU for a 
supersonic operator in Eq. (11) is different from 
the definition for the subsonic operator of Eq. (12). 

The evaluation of the subsonic operator in 
Eq. (12) requires velocity potential ( <t>) values at 
i+l and i+2 planes from the previous (n) relaxation 
cycle to compute the density. The Initial and 
Boundary Conditions section of this paper pre- 
scribes a method to start the first relaxation 
cycle of the subsonic bubble calculation. 

B. Treatment of |p in Eq. (1) 


When the point is elliptic , the density biasing 
is defined by 



where v = maxlO, 1 -— ). As before, the super- 
\ q 2 / 

script n+1 denotes the current relaxation cycle for 
a subsonic bubble calculation. Note the difference 

in the definition of v and u. The density biasing 
in the cross flow direction n is turned off when 
the total velocity q is less than the speed of 
sound a, just as in the marching c direction 
( Eq . (13)). The implicit treatment of V in the 
marching subsonic operator of Eq. (15) is the same 
as that of the supersonic part, explained in 
Ref. 5. 


Referring to the Fig. 3a molecule, 



supersonic 


where 

e i+l 


+ (1 


- 9,. 


i+l 



marching 

subsonic 


if 



Hi) < o 

i+l 


(supersonic 

point) 


A similar procedure is implemented for the 
(pjj term in Eq. (1). 

D. Implicit Factorization Algorithm 

Combining the various terms of Eq. (1) as repre- 
sented by Eqs. (10)— (17) together with the terms 

arising from |p jj will result in a fully implicit 

model. This is solved using an approximate factor- 
ization implicit scheme. After some rearrangement 
of the terms, the factored implicit scheme 
becomes 
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A$ = R . 


08) 


. 1 _3_ ^22 _3_1 
F 3n T" 3n| 

The density p appearing in Eq. (18) can be either 
H or B depending on the sign of ~~j as 

illustrated in Eq. (15) . 

Equation (18) has the form 

L c L n (A<(.) = R (19) 

and it is implemented as follows: 

L 5 (A4>)* = R L (A<(.) = (Ab)* * i+1 = <|>. + A$ . (20) 

The various quantities appearing in Eq. (18) are 
given by 

A > ' & (*" ' $ 


Inside an MSR, as in Fig. 2, when Eq. (12) is 
applied at an (i+1) grid point, information on q. 


i+2 


is required to form the density p and various 
derivative terms. For the first relaxation pass, 
an initial estimate for quantities in the (i+2) 
plane is prescribed in the following manner: 




MSR 

operator 


needs 

initial 

estimate 


' (j (a l l^c + a 12^n + a l 3 (t 5 ) )j i + 1 


( 22 ) 


In Eq. (22), sonic conditions are assumed at 
(i+2) for the first relaxation pass: 


i+2 


u i+2 = q * (a n ) i+2 


(23) 


The sonic values p* and q* are purely a func- 
tion of the freestream Mach number M^. Also, P i+ -| 

in Eq. (22) is initialized to be P i . 

For the second relaxation cycle and onwards 
(n > 1), the conditions from the previous relaxa- 
tion cycle are used: 


A 3 = e i 


r p i 


+ 

CD 

1 

1 

(H, 


■(•Cl 

i 

[ (24) 

r 




n+1 . 

n 

1 

p i 

h+i 

(•»-“) 

-°- 6 i+i ) 

pa 13\ 

J /i+1 

p i+l " p i+l ) 

Boundary Conditions 



Ac 0 = ^i+2 ‘ c i+l 


* = ? i+l - ?i 


( 21 ) 


At a solid boundary, the contravariant velocity 
V is set to zero. Exact implementation of V = 0 in 
the implicit treatment of Eq. (18) is described in 
Ref. 4. 


and the right-hand side term R consists of various 
known quantities. 

If the flow field does not contain an embedded 
MSR or TSR, the implicit factored algorithm of 
Eq. (18) performs a pure marching procedure start- 
ing from an initial known data plane. In this 
situation, there is no need to go back to the 
upstream starting plane and iterate the solution. 
However, If a subsonic bubble is present (between 
planes AB and CD in Fig. 2), then the solution 
procedure of Eq. (18) performs a relaxation method, 
and iterates for the elliptic subsonic bubble to 
converge (superscript n in Eqs. (12), (13), and 
(17) refers to the relaxation cycle counter). 

E. Initial and Boundary Conditions 

Initial Conditions 

For a pure supersonic flow, initial conditions 
need to be prescribed only at the starting plane. 
Usually, the starting plane is set close to the 
apex of the configuration solved for, and conical 
solutions are prescribed. 


The outer boundary is set away from the bow 
shock and the freestream velocity potential is 
imposed along that boundary. All discontinuities 
in the flow field are captured. The precise density 
biasing activator v, based on the characteristic 
theory, allows for sharp capturing of shocks in 
the flow. 

Behind the trailing edge of a wing, a wake 
model is imposed. Figure 4 shows a schematic of a 
wake model. At a point P lying on the wake, the 
boundary condition is that there is no jump in 
the pressure across the wake, i.e., (p p -Pg) = 0. 

In the full potential (isentropic) formulation, 
this translates into the condition that the jump 
in density (p p -Pg) is zero, or the jump in the 

total velocity q is zero ((q p -qg) * 0). The jump 

in q across the wake is set to zero in an approxi- 
mate manner in the following way. 

First, compute the jump in the potential $ at 
the trailing edge point P 1 , and maintain that jump 
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to,?)p ■ G ^rjn^p ~ 0 


Fig. 4. Wake boundary condition. 


constant along the line P'P in Fig. 4. At the wake 
point P, Eq. (1) is not valid. Instead of solving 
Eq. (1), 4> nn = 0 is satisfied at the wake point P, 

to achieve the condition (<t> ) - U ) = 0. Incor- 

' n 'p ' n'q 

porating a constant jump in 4> along P'P ensures 



is approximately set to zero, yielding the necessary 
wake boundary condition. The Results section in 
this paper presents a calculation performed for a 
realistic wing-body-wake fighter model and shows 
an excellent matching of the pressures across the 
wake, using the above wake boundary condition. 


F. Grid System 

The transformation from the physical space 
(x,y,z) to a body-fitted computational space 
(c.n.E) is performed numerically at each constant 
t plane by using the elliptic grid generation 
technique of Ref. 6. Once the grid is generated, 
all the metric terms a.., in Eq. (2) and the 

• J 

Jacobian J in Eq. (3) are computed by numerical 
differentiation. As described in Ref. 5, a free- 
stream error subtraction is performed at each grid 
point to account for any improper metric 
cancellation. 

G. Density Biasing Summary 


This section summarizes in a tabular form the 
type-dependent density biasing procedures incorpor- 
ated in this paper to generate the proper artificial 
viscosity. 
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Term 

Total 

Supersonic 

( a n-p-) <0 * 

Marching 

Subsonic 

Total 

Subsonic 

( a irp-) >0 » 


q>a 

q>a 

q<a 

pU in 

t-direction 

Upwind 

differ- 

encing 

Eq. (11) 

Dens i ty 
biasing 
based on 

{'-$) 

in 

Eq. (13) 

Shut off 

density 

biasing 

P V, pW 
in n,E 
directions 

Density 
biasing 
based on 

I 1 ‘ a 22 Jr) 

Density 
biasing 
based on 

•("?) 

Shut off 

density 

biasing 


i 1 * a 33 Jr) 

p in 

Eq. (17) 



P in 
Eq. (16) 




Table 1. Sumary of type-dependent density 
biasing procedure 


IV. Results 

As illustrated in Fig. 2, supersonic marching 
calculations are performed from the nose until an 
embedded MSR forms. In Fig. 2, the plane AB is the 
last supersonic marching plane preceding the sub- 
sonic bubble and forms the upstream computational 
boundary for the relaxation calculation. For the 
first relaxation pass through the subsonic bubble 
region, 0 i+1 in Eq. (10) is set equal to and 

( p U)i + 2 = P*q*. From the second relaxation cycle 
on, 0 i+ -|. 6,j , and (pU) i+2 are computed according to 

their definitions. A typical supersonic flow with 
a subsonic bubble calculation required at most 
only four relaxation cycles (iterating back and 
forth between planes AB and CD) to obtain a 
converged location for the bubble. The initial 
guess, based on the sonic conditions p*q*, worked 
out very well for all the subsonic bubble cases 
presented in this paper. The (n.E) marching plane 
can be any arbitrary surface, but for convenience 
was chosen to be a constant x-plane. 

The step size in the marching direction, c, 

for the supersonic part I " ~t) < 0, q > a) , 

was automatically chosen by setting the Courant 
number^ to be around 5. Once the MSR forms, the 
eigenvalues become complex, and the step size can 
not be computed based on a specified Courant num- 
ber. For marching planes containing the MSR/TSR, 
the step size was specified into the code depending 
on the geometry variation. When geometry changes 
were drastic (region of emergence of a wing from 
a fuselage), usually a smaller step size At was 
required (as small as .003^.005 for a total 
length of one) to properly account for rapid 
changes in the flow. Once the MSR/TSR is fully 
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captured and the flow becomes supersonic again, the 
step size selection once again becomes based on the 
Courant number. For a pure supersonic flow all the 
way, the entire calculation could be performed 
using 40 planes or less (Ac > .025). However, once 
an MSR or TSR is present, the total number of 
C planes in the calculation could go as high as 150. 

Figure 5 shows the surface pressure distribu- 
tion in the axial direction on the upper (6 = 0, 
lee side) and lower (6=180°, windward side) plane 
of symmetry for a developed cross-section forebody 
geometry reported in Ref. 7. At 1.7 and 
a = -5“ , the lee side has an embedded MSR which 
required use of the relaxation operator in Eq. (10). 
A pure supersonic X-marching for this case would 
have failed without the MSR treatment described in 
this paper. 

Figure 6 shows the axial pressure distribution 
on the same geometry for M x = 2.5, a = -5°. This 
case has no embedded subsonic bubble and the entire 
flow remains supersonic in the marching X-direction. 
The same code is used to generate the results of 
Figs. 5 and 6. The code automatically turns on the 
subsonic bubble relaxation treatment without the 
user having to play a role. 

Figure 7 shows the circumferential pressure 
distribution for the same developed cross-section 
forebody at M^'1.70, a=-5°, and x/8. = 0.28. The 
embedded MSR thickness is the largest at this 
axial station. The extent of the subsonic bubble 



Fig. 5. Axial surface pressure distribution for 
a developed cross-section forebody; 
H.-1.7, a=-5". 
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Fig. 6. Axial surface pressure distribution for 
a developed cross-section forebody; 
M 00 =2.5, a = -5° . (Fully supersonic 
flow in the marching direction.) 

is marked in Fig. 7. The results of Fig. 5 and 
Fig. 7 only exhibit MSR, and TSR is not present. 

To simulate both the MSR and TSR, a case was 
constructed with a triple cone geometry. At 
M = 1.6 and a = 0° , this case exhibited the presence 
of both the MSR and TSR, as shown in Fig. 8. The 
density biasing is activated in MSR, and turned off 
automatically in TSR by Eqs. (13) and (17). 

For the case of Fig. 9, the pressure distribu- 
tion in the n direction (away from the body) at an 
axial station containing the MSR and TSR is shown. 
The sonic C* at this Mach number is 0.695. Outside 

of the bow shock, the pressure coefficient C goes 
to zero as shown in Fig. 9. " 

The axial surface pressure and local Mach num- 
ber distribution for the triple cone geometry of 
Fig. 8 are shown in Fig. 10. The MSR and TSR 
locations are marked. The local Mach number goes 
below 1 inside the TSR. 

Figure 11 shows a supersonic fighter configura- 
tion with a wing sweep of around 48°. At a free- 
stream Mach number of 1.6 and a = 5°, the leading 
edge of the wing exhibits an MSR. To solve the 
flow field over such a fighter configuration, one 
needs to use the embedded subsonic bubble treat- 
ment. Figures 12 and 13 show the surface pressure 
at various axial stations along with respective 
grid distribution for the wing-body geometry. The 
leading edge behaves like a subsonic leading edge, 
producing a suction peak. For this case, the MSR 
starts around x = 0.4. Figure 14 shows the pressure 
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Fig, 7. Circumferential pressure distribution 
for a developed cross-section forebody; 
M^ = 1.7, a= -5°, x/i = 0.28. 
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Fig. 8. Marching subsonic and total subsonic 
regions embedded in a supersonic flow 
over a typical fuselage-canopy 
geometry; M^'1.6, a=0 6 . 

distribution for the fighter configuration of 
Fig. 11 at an axial station x/2. = 0.85, where a 
wake sheet is present. The grid distribution goes 
around the wake sheet just like a wing-body case. 
The approximate wake model described in the paper 
seems to provide the correct zero pressure jump 
condition across the wake, as seen in Fig. 14. 
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(v - y b > 


Fig. 9. Pressure distribution away from the 
fuselage-canopy geometry in the 
y-direction; = 1.6, a=0°. 



Figure 15 shows the pressure distribution Fig. io. Axial surface pressure and Mach 
right along the leading edge of the fighter con- number distribution on a typical 
figuration. At the axial location where the wing fuselage-canopy geometry. 
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Fig. 11. Supersonic fighter with an embedded 
marching subsonic region near the 
leading edge. 

emerges from the body, the C rises sharply, 
creating an MSR. p 

For the wi ng-body-wake model, at ^=1.6, 
a = 5° , the overall lift coefficient (without 

the tail) came out to be 0.298, while the tunnel 
data was 0.277. The overall drag coefficient C D 
from the code was 0.0462 (after accounting for u 
a friction drag coefficient value of 0.01077), 
compared to the tunnel data of 0.0457. 

Figure 16 shows the drag prediction capability 
of the full potential code by demonstrating it on 
a double wedge delta wing at M oo =1.62. At this 
Mach number, the leading edge exhibited the 
presence of an MSR for sweep angles less than 60°. 
A pure supersonic marching code would not have 
worked for this case. The drag calculation from 
the full potential code compared very well with 
the experimental data available in the Princeton 
Series. 



Fig. 12. Pressure distribution on a fighter- 
like configuration, M^l.6, a =5°. 


V. Conclusions 

A nonlinear full potential method has been 
developed to treat supersonic flows with embedded 
subsonic regions. A conservative switching scheme 
is employed to transition from the supersonic 
marching algorithm to a subsonic relaxation 
procedure. The theory of characteristic signal 
propagation plays a key role in activating various 
density biasing procedures to produce the neces- 
sary artificial viscosity. The method has been 
shown to produce results which were hitherto not 
possible using a pure supersonic marching scheme. 
Currently, work is progressing to treat complex 
configurations with vertical tails, like the one 
shown in Fig. 17. The concept of density biasing 


will also be modified in the future to a flux 
biasing procedure described in the Appendix. 


Appendix 

Flux Biasing Procedure 

O 

Based on the work of Hafez and Osher , it is 
possible to modify the density biasing concept to 
a flux biasing procedure. 

Consider the term — |p jj in Eq. (15). The 
density biasing procedure defines p to be 
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Fig. 13. Pressure distribution on a fighter- 
like configuration at x = 0.55, 

M^l.6, a= 5°. The leading edge 
is in a marching subsonic region. 
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Fig. 14. Grid and pressure distribution in 
the wake region of a fighter-like 
configuration, M =1.6, a =5°, 
x = 0.85. 
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Fig. 16. 


Drag prediction for a double wedge 
delta wing at fQ= 1 .62, a = 0°. 

For sweep angles less than 60°, 
the leading edge has a marching 
subsonic flow. C D = C n 
u min u f 



nviscid 


Fig. 15. Pressure distribution along the lead- 
ing edge of a fighter configuration 
at fQ'1.6, a=5°. The leading edge 
has a marching subsonic region (MSR). 
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In the flux biasing technique, it will be 
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JL 

3n 


Jq”)Pq - Ari 


$ n (pq) 


I (A2) 

Uj+H 


where 


(pq)_ * 0 

= (pq) - p*q* 


if 

q < a 

if 

q > a 


7. Townsend, J.C., Howell, D.T., Collins, I.K., 
and Hayes, C., "Surface Pressure Data on 

a Series of Analytic Forebodies at Mach 
Numbers from 1.7 to 4.50 and Combined 
Angles of Attack and Sideslip," NASA TM 
80062, June, 1979. 

8. Hafez, M., Whitlow, W. , Jr., and Osher, S., 
"Improved Finite Difference Schemes for 
Transonic Potential Calculations," 

in preparation. 


where p* and q* represent sonic conditions, a is 
the local speed of sound, and (pq) is the flux. 
When the flow Is purely subsonic, the flux biasing 
is turned off automatically. 
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Abstract 

A two-dimensional elliptic grid solver Is presented 
and its application to various three-dimensional con- 
figurations, both internal and external, is demon- 
strated. The method uses proper forcing terms to clus- 
ter grid points near boundaries with a specified grid 
spacing and allows grid lines to intersect the bound- 
aries at a specified angle. By segmenting the region, 
grid results are generated for sharp leading edged con- 
figurations and wing-vertical tail combinations. 


Introduction 

An integral part of the Computational Fluid 
Dynamics work is the development of numerical grid gen- 
eration procedures for a body-fitted coordinate system 
as a practical way to apply boundary conditions. Var- 
ious schemes are available to achieve such a body- 
fitted system, including conformal mapping, 1 algebraic 
schemes, 2 elliptic grid sol vers, 3 * s and the hyperbolic 
solver. 6 The method commonly in use is the elliptic 
grid solver which solves a set of Poisson equations 
with appropriate right hand side forcing terms to 
achieve two main features: 1) to cluster points opti- 
mally near the boundary and in regions of high gradi- 
ents in flow, and 2) to force grid lines to intersect 
the body surface and other computational boundaries in 
a nearly orthogonal fashion. In elliptic grid solvers 
the quality of the grid distribution critically depends 
on the choice of the forcing terms. 

The purpose of this paper Is to illustrate the 
effective use of one such elliptic solver with suitable 
forcing terms, in developing computational grid for a 
wide variety of -oth internal and external three- 
dimensional problems. Examples to be shown in this 
paper Include grid generation for a three-dimensional 
inlet (both interior and exterior of the inlet), exter- 
nal grid for a sharp leading edged wing-body and wing- 
vertical tail combinations, and grid for multiple- 
connected regions of an arrow wing-body configura- 
tion. The key result in the paper is the grid genera- 
tion for a three-dimensional inlet by using just a two- 
dimensional grid solver. 


The elliptic grid solver to be presented here Is 
somewhat similar to that of Ref. 5, but differs consid 
erably in the implementation of a proper procedure for 
the evaluation of the forcing terms. 


Formulation 

Thompson et al. 3 have proposed the following inho- 
mogeneous Laplace equations as the governing system for 
grid generation 


? xx + «y> ‘ PU>n) 

n xx + n yy » Q(E,n) . (1) 

Equation (1) represents the transformation from the 
(x,y) physical domain to the (E.n) computational 
plane. The right hand side of Eq. (1) defines the 
forcing terms to achieve a desired grid distribution. 
Various grid generation methods 3 * 5 differ from each 
other In the choice and implementation of appropriate P 
and Q source terms to provide desired grid proper- 
ties. The purpose of this paper is to describe an 
efficient and automated procedure for choosing the 
proper forms of P and Q to achieve grid clustering and 
grid orthogonality near surfaces. 

Interchanging the role of dependent (c,n) and inde- 
pendent (x,y) variables, In the transformed space, 

Eq. (1) takes the form 

a ' 2sx ?n + TX nn * + Qx n } 

“ y « ' 2By 5n + Ty nn = * j2(Py E + Qy n ) (2) 

where a = x 2 + y 2 , 5 = + y^, y = x| + y| and 


Choosing the following forms for P and Q 
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P = J 2 4(E,n)a 

Q » J 2 4(E.nh , {3} 


Equation (2) takes the form 


• C X EE + * X E } ' 2iJx En + Y(x nn * * 0 

a(y a ♦ « e ) - 2sy 5n -r(y„ n ♦ w,) « 0 . (4) 

Equation (4) contains two free parameters, 4 a^yJ< 4 . 
Referring to Fig. 1, the role of these parameters is to 
enforce two conditions that the grid lines have^to 
satisfy. 1) Maintain a given aS grid spacing in the 
transverse direction (n) at the surface n = n b , and 2 ) 
provide a specified intersecting angle 6 between the 
surface (n = nJ and the transverse coordinate (e * 
constant). Solving for 4 and 4 from Eq. (4), one gets 


( *n ) n«n b " S n ( - X E cos 6 * h s1n e )/^p~y*" 

( 9 ) 

k„)„«n b * S n ( ‘ y E cos 6 + X E s1n •) / ^f T 

where e and are specified tjy the user. Once the 

grid points along n = nb (body surface) are prescribed, 
the differential operators x^, y^, x^, y^ anoearing 

In Eq. (S) car be evaluated at the body surface 
(n * n^) using the following relationships (refer to 
Fig. 1 for notation) 

V, " ( * 7x l + 8x 2 ‘ *3)^ 1,2 - 3 (* J ,) t ,« Tlb / 4 i 

( 10 ) 

y nn “ *' 7y l + 8 *2 * * 3 Vn=n b /AT1 


♦ * ’7 <VtC * VeE 5 + S (x En y n * y En X n> 

o7 ^n x nn ~ X n y nn^ 

♦ " 1 (x nn y E ' y nn X E } “ K ( We * y En X E } 

+ ^*EE y E * y EE*E ) • 

Equation (5) is valid at every grid point. However, 
the evaluation of 4 and 4 requires values for all the 
derivatives V * ? , y„, y ? . x^, y^, y^, x^ 

and y . One can, however, solve for 4 and 4 on the 

surface (n - n b ) in Fig. 1, by prescribing the aS and e 
values in the following manner: 

1. Specify the Intersecting angle 8 between the trans- 
verse grid line (e * constant) and the body surface 
(n = n b ). This condition can be written as 


VE-Vn « |ve | |vn| cos e 


( 6 ) 


In the transformed space, Eq. ( 6 ) takes the form 


where (*„)„.,, and (y ) are obtained from Eq. (9). 
'b b 


Using the derivative values from Eqs. (9) and (10), the 
expressions for 4 and 4 given by Eq. (5) can be pre- 
cisely computed along the boundary n ■ ok. The same 
procedure Is repeated along all the boundaries of the 
computational domain to obtain the boundary variation 
of 4 and 4 . Once 4 and 4 values are known along the 
computational boundaries, their values at interior mesh 
points are computed by interpolation. Usually 8»n/2 is 
used in Eq. (7) to achieve nearly orthogonal grid lines 
at the boundaries. The use of transverse direction 


derivatives x_ 


V Vi- 


no 


x . and y _ in the 
nE nE 


evaluation of 4 and 4 In Eq. (5) makes the procedure of 
this paper considerably different from the technique of 
Ref. 5. The 4 and 4 values along boundaries are con- 
tinuously updated after each grid relaxation cycle, in 
a manner similar to the procedure in Ref. 4. 


The procedure described so far generates a two- 
dimensional grid in the x,y plane. If the boundary, 
say n « n b in Fig. 1, is not confined to a two-dimen- 
sional plane, that is, along n « n b , all x, y and z are 
specified as shown in Fig. 2. Then, one can generate a 
warped computational plane that contains the bounoary 
n ■ n b in the following manner: 

1) First solve for x,y from Eq. (4). This pro- 
vides x,y grid in the projected plane, shown 
in Fig. 2. 


X E X n + y E y n “ * ^ X n + y n^ X E + cos 8 . 

( 7 ) 


2) Then solve the Laplace equation 


z xx + z yy 


0 


( 11 ) 


2. Specify the grid spacing &S at the boundary 
(n = n b ) in the transverse direction (e = constant). 
This condition can be stated in the form 


S 


n 


= 7x2 + H • 
n n 


( 8 ) 


Using Eqs. (7) and (8) and the relationship 

VE * Vn = |7E||^n|sin 6, one can write the following 

relationships along the body surface n = n b . 


with z prescribed along boundaries, 
transformed space, Eq. (11) -ecomes 

“< Z EE + * Z E } ' 28 z En- y(z nn + * Z n ) 


In the 

= 0 . 
( 12 ) 


Solution to Step (1) above provides values of a, 6,^,4 
and 4 that are required in solving Eq. (12). 

The grid results to be presented in this paper are 
obtained by solving Eqs. (4) and (12) simultaneously. 
The method has been effectively used to generate both 
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the Internal and external grid for a three-dimensional 
inlet system, and also the external grid for varieties 
of sharp leading edged wing-body combinations. A typi- 
cal grid calculation required approximately 100 relaxa- 
tion cycles to converge the residual to 10' 8 level. 

Results 

First external grid results are presented for sharp 
leading edged configurations, wing-vertical tail com- 
binations and then grid results for a three-dimensional 
inlet are shown. 

External Grid 

In the present study, the external grid for a wing- 
body combination is generated for use in a supersonic 
flow calculation by a marching technique. 7 This re- 
quires grid generation in every marching plane, which 
can be either a constant z plane or a spherical plane. 
Figure 3 shows the schematic of a sharp leading edged 
wing-body cross section. To be able to accommodate 
sharp edges in the grid generation routine, a cut A-B 
in Fig. 3 is made and the regions 1-2-8-A-l and 
3-A-B-4-3 are treated separately with a given grid 
distribution along the cut A-B. 

Figure 4 shows a double-wedge delta wing and the 
grid generated at various axial cuts of the geometry. 
Figure 5 shows the pressure contours at a typical axial 
station of the geometry in Fig. 4 indicating the pres- 
ence of a shock and expansion fans at M = 1.92, a « 

0°. Figure 6 shows the minimum drag results for zero 
lift and comparison with experimental data. The suc- 
cess of a good flow calculation such as the one shown 
in Fig. 6 critically depends on the quality of the 
grid. 

Figure 7 shows the grid generated for an arrow 
wing-body combination. The cross section starts 
initially as a wing alone, then transitions to a wing- 
body, and finally becomes a detached wing-body. The 
grid generation routine is general enough to accommo- 
date such a shape change in the axial direction. Flow 
calculations obtained for the arrow wing-body combina- 
tion using the grid of Fig. 7 are shown in Fig. B. 
Improvements to the wake treatment for the methodology 
of Ref. 7 are currently underway. 

Figure 9 shows the schematic of a typical fighter 
configuration and grid results at various axial loca- 
tions. By segmenting the entire region, the present 
grid solver can be effectively used to generate grid 
around wing-vertical tail combinations. Figure 10 
shows the gridding for a sharp leading edged Sears- 
Haack wing-body combination. 

Three-Dimensional Inlet Grid 

Figure 11 shows a schematic of a three-dimensional 
inlet embedded in a global computational domain. The 
objective is to develop a body-fftied coordinate system 
both on the exterior and the interior of the inlet. 

The le-ding edge of the inlet (1-2-3-4-5-6-7-1 in 
Fig. U) is highly swept and curved. The intent is to 
perform a Navier-Stokes calculation for that geometry, 
which will require a clustered grid near the leading 
edge as well as near the walls of the inlet. To create 
a clustered body-fitted as well as leading edge fitted 
grid system, it was decided to gradually rotate the 
upstream constant z plane such that there will be 
leading edge plane containing the entire leading edge 
1-2-3-4-5-6-7-1. This warped leading edge plane would 


then be gradually rotated to a constant z downstream 
plane. In order to achieve this, first series of loops 
were created starting from the leading edge to the 
downstream plane as shown in Fig. 12. For each loop, 
the inner and outer inlet wall geometry in terms of 
(x.y.z) was prescribed as shown in Fig. 13. Then, for 
each loop a warped plane containing the body-fitted 
grid system was generated using the x,y grid solver 
procedure of Eq. 4 and then a;“z* value at each grid 
point from Eq. (12). 

Figure 14 shows a loop near the leading edge of the 
inlet, a perspective view of the warped body-fitted 
grid surface, and a constant z projected plane showing 
the details of the interior and the exterior grid. 

Figure 15 shows an enlarged view of the upper and lower 
left hand corner region of the inlet. Figure 16 shows 
a loop at a mid-station location of the inlet clearly 
indicating the thickness of the inlet wall. Also shown 
in Fig. 16 are the perspective view of the warped body- 
fitted surface and a z-plane projected grid. Figure 17 
shows the grid at a final downstream plane for the inlet. 

Figures 11 to 17 represent three-dimensional inlet 
grid work obtained purely by using two-dimensional grid 
solvers based on Eqs. (4) and (12). Currently, this 
inlet grid is being incorporated into a three- 
dimensional Navier-Stokes code. 

Conclusions 

This paper presents an elliptic grid solver in two 
dimensions with proper forcing terms to obtain certain 
desired grid qualities such as grid clustering and 
orthogonality near body surfaces. The procedure has 
been successfully applied to various three-dimensional 
external Sharp leading edged configurations as well as 
for a three-dimensional inlet. 
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Fig. 1 Notation for the grid solver. 



PLANE (a) 
PROJECTION 
ON U, y) PLANE 
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Fig. 3 Schematic of a sharp leading edged wing-body 
cross section. 


Frg. 2 Grid generation on a curved surface. 
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Gridding for a double wedged delta wing. 
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Fig. 5 Pressure contours at x/t ■ 0.4 for the delta 
wing at = 1.92 .and zero degree angle of 
attack. 


M. . 1.92 



(TAN E/TAN «) 

Fig. 6 Minimum drag coefficients for the delta wing 
and comparison with experimental data. 
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Fig. 7 Gridding for an arrow wing-body combination. 



Fig. 8 Flow calculations for the arrow wing-body 
configuration at * 2.96 and a = 10.01°. 
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Fig. 9 Gridding for a typical fighter configuration 
with vertical tails. 
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Fig. 14 Body-fitted grid system on a warped plane near 
the inlet leading edge. 



b) LOWER RIGHT HAND SIDE CORNER 


Fig. 15 Enlarged views of the body-fitted coordinate 
system near corners of the inlet. 
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Fig. 17 Body-fitted grid system at the final 
downstream plane. 


Fig. 16 Body-fitted grid system at the mid station of 
the inlet. 
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Fig. 13 Inner and outer wall shapes of the inlet at 
various loop locations. 


OUTER WALL 
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